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Sensitivities  and  Variances  for  Fitted  Parameters  of  Spheres1 * 3 

Christoph  WitzgalT 
Marek  Franaszek^ 


Abstract 

Experimental  data  have  been  gathered  by  applying  3D  imaging  systems,  such  as 
LID AR/LAD AR( -LIght/LAser  Detection  And  Ranging)  instruments,  to  spherical  objects.  This 
report  provides  a compilation  of  the  statistical  and  analytical  procedures  to  be  used  for  an 
analysis  of  these  data  to  be  reported  separately.  This  analysis  will  be  based  on  two  different 
nonlinear  least-squares  approaches  to  modeling  objects,  directional  and  orthogonal  fitting.  It  is 
proposed  to  estimate  variances  of  fitted  sphere  parameters  directly  from  their  sensitivities  to  data 
perturbations  rather  than  to  follow  the  common  prescription  of  linearization.  The  sensitivities  are 
determined  by  implicit  differentiation  of  error  gradients.  Detailed  descriptions  of  the  directional 
and  orthogonal  fitting  methods  are  set  forth  as  applied  to  spheres  in  a scanning  environment.  In 
particular,  the  report  furnishes  closed-form  expressions  for  those  derivatives  of  the  respective 
error  functions  which  are  needed  for  the  calculation  of  the  parameter  sensitivities  with  respect  to 
the  full  set  of  control  variables. 

Key  Words 

3D  imaging  systems,  directional  fitting,  least  squares,  LIDAR/LADAR,  orthogonal  distance 
regression,  parameter  variance,  parameter  sensitivity,  sphere,  statistical  modeling 


1 Introduction 

A frequent  task  is  to  determine  the  shape  characteristics,  size,  location  and  pose  of  physical 
objects  for  purposes  of  identification,  location,  registration,  and  calibration  of  coordinate  frames. 
Such  tasks  are  needed,  among  others,  for  quality  control  in  manufacturing,  determination  of  “as- 
built”  structures,  construction  automation  and  site  monitoring,  e.g.,  [1,2], 

A common  approach  to  these  tasks  is  to  acquire  3D  coordinates  of  data  points  considered  to  lie 
on  the  surface  of  targeted  objects.  3D  Imaging  Systems,  which  include  “line-of-sight" 
LID  AR/LAD  AR(  ^LIght/LAser  Detection  And  Ranging)  devices,  are  increasingly  used  for  this 
purpose.  The  latter  instruments,  in  particular,  are  capable  of  fast  generation  of  large  amounts  of 
data  points  or  “point  clouds”.  They  scan  an  object  by  emitting  laser  pulses  and  processing  return 
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signals  in  order  to  determine  the  distance  traveled  and  thus  determine  the  distance  or  “range” 
between  the  instrument  and  the  point  of  impact  — presumably  - on  the  object.  The  device  keeps 
track  of  the  “bearings”  such  as  azimuth  and  elevation  angle  at  which  each  particular  signal  was 
emitted.  This  process  of  data  acquisition  suggested  the  use  in  this  work  of  polar/spherical 
(“angle-angle-range”)  coordinate  systems  for  representing  data  points.  Also,  transformation  to 
Cartesian  coordinates  would  introduce  correlation. 

Once  the  point  cloud  corresponding  to  an  object  has  been  determined,  a computational  process  is 
required  to  extract  the  desired  features  of  the  object  from  this  data  set.  In  typical  applications,  a 
mathematical  “model"  is  specified,  based  on  features  characteristic  for  a class  of  objects.  The 
model  is  “parameterized”,  that  is,  it  is  defined  with  the  help  of  parameters  that  determine  these 
characteristics.  Choosing  values  for  these  parameters  will  result  in  the  mathematical  description 
of  a surface  to  represent  a “virtual  object”,  which  may  then  be  compared  to  an  image  of  the  real 
object  as  provided  by  the  point  cloud.  By  adjusting  the  model  parameters  so  that  the  virtual 
object  moves  into  a location  that  optimizes  its  proximity  to  the  object  image,  the  desired 
characteristics  such  as  location,  pose,  size  or  shape  of  the  object  are  found  within  the  coordinate 
frame  of  the  point  cloud.  This  permits  determining  the  geometric  relationship  between  that 
object  and  other  objects  or  features  that  are  also  represented  in  the  point  cloud  frame.  If  this 
frame  registers  to  an  established  ground-truth  frame,  then  absolute  measurements  of  location, 
pose  and  shape  can  be  extracted. 

Approaches  to  object  modeling  may  employ  the  powerful  “Iterative  Closest  Point  (ICP)”  method 
[3],  or  the  “Hough  Transform”,  e.g.  [4],  Present  work  focuses  on  the  extensively  used  “Fitting” 
paradigm,  which  is  based  on  minimizing  a specified  error  function  or  on  maximizing  likelihood. 
The  reader  may  want  to  consult  texts  on  “Statistical  Models”  such  as  [5-7]. 

Of  particular  interest  are  two  least-squares  based  approaches,  “orthogonal”  and  “directional” 
fitting.  Orthogonal  fitting,  also  referred  to  as  “Orthogonal  Distance  Regression  (ODR)”  [8,  9],  or 
“Geometric  Fitting”  [10],  is  a commonly  used  and  widely  commercialized  method.  Many 
publications  [10-19]  discuss  its  application  to  the  fitting  of  spheres  or  circles.  The  alternate 
approach,  “directional  fitting”,  has  been  proposed  and  discussed  [20,  21]  for  data  acquired  by 
scanning  from  a single  instrument  position.  Here,  the  orthogonal  (closest  Euclidean)  distance  to 
the  virtual  object  has  been  replaced  by  the  distance  in  the  direction  of  the  scan  by  which  the  data 
point  had  been  acquired.  While  computational  concerns  dominate  much  of  the  field,  our  interest 
here  is  in  statistical  and  metrological  issues. 

At  the  core  of  this  report  is  an  approach  to  determining  the  “sensitivities”  of  fitted  model 
parameters,  in  general,  and  for  spherical  models,  in  particular.  The  report  is  also  preparatory  to 
an  experimental  application  of  different  fitting  methods  and  their  statistical  evaluation  [22]. 
Given  specified  variances  for  range  measurements,  the  derivation  of  variances  for  fitted  sphere 
centers  is  a major  interest.  In  Chapter  2,  the  general  fitting  paradigm,  based  on  the  concept  of  an 
error  function,  is  described  along  with  the  general  computational  formalism  for  calculating 
parameter  sensitivities.  These  sensitivities  will  be  used  to  estimate  parameter  variances.  In 
dealing  with  spherical  models,  this  approach  is  of  necessity  more  general  than  the  common 
nonlinear  least  squares  approach  based  on  linearization  and  homoscedacity.  Chapters  3 and  4 are 
dedicated,  respectively,  to  orthogonal  and  directional  fitting  of  spheres  in  a scanning 
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environment.  Closed  forms  of  the  derivatives,  needed  for  calculating  sensitivities,  are  reported. 
The  Appendix  will  feature  detailed  derivations  of  the  reported  formulas  so  as  to  enable 
verification. 


2 The  Fitting  Paradigm 
2.1  Error  Function 

Once  a parameterized  model  has  been  selected,  it  is  natural  to  ask  for  parameters  that  minimize 
the  extent  to  which  the  point  cloud  deviates  from  the  resulting  virtual  object.  The  hope  is  that 
such  an,  at  least  locally,  optimal  virtual  object  provides,  within  the  coordinate  frame  of  the  data 
points,  an  accurate  representation  of  the  actual  object.  Fitting  a 3D  model  of  a sphere  of  a 
Cartesian  center  C = [X,T,  Z]  and  known  radius  R may  be  accomplished  by  specifying  an 
“error  function” 


(2.1.1) 


to  be  minimized  by  varying  the  model  parameters  X,Y,  Z,  while  the  data  variables 


dr(pr6r  are  coordinates  of  points  P = \dt  q>t  0t]  to  be  measured.  The  following 


discussions,  however,  should  not  be  construed  as  pertaining  only  to  this  special  scenario,  but 
rather  as  representative  of  full  generality.  In  particular,  the  data  may  also  be  Cartesian,  or  not  be 
coordinates,  at  all.  The  choice  of  the  error  function  should  be  such  that  it  produces  only 
nonnegative  values.  A minimum  of  zero  should  indicate  a perfect  fit.  An  error  function  E thus 
furnishes  a model  description. 

Given  an  actual  data  set  of  measurements 


the  parameter  values 


= A 


are  thus  determined  by  minimizing  the  expression 


for  the  variables  X ,Y,  Z , given  the  coordinate  values  of  the  data  points  P " . 


A common  approach  to  constructing  error  functions  is  to  assign  an  individual  error 
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to  each  stipulated  point  P = [di  (pi  0)  ],  and  to  minimize  the  sum  of  squares 


(2.1.2)  E = J>,2 . 

1=1 

The  orthogonal  and  directional  fitting  methods  are  based  on  this  Nonlinear  Least  Squares  (NLS) 
concept  [17,  19].  In  both  methods,  each  point  P is  assigned  a “theoretical  point”  or  “model 

point”  P located  on  the  proposed  virtual  object.  That  theoretical  point  is  seen  as  the  desired 

“correct”  point,  and  the  Euclidean  distance  between  the  two  points  is  considered  the  individual 
error 


e.  - 


P -P 


of  the  data  point  with  respect  to  the  current  location  and  shape  specification  of  the  virtual  object. 

In  orthogonal  fitting,  the  theoretical  point  P is  chosen  as  a point  that  lies  on  the  virtual  object 
and  is  closest  to  the  data  point  P in  terms  of  Euclidean  distance.  In  the  3D  imaging 
environment,  however,  the  data  point  P is  considered  to  lie  on  a particular  “scan  ray”  or  “line- 
of-sight”,  which  emanates  from  the  instrument  position. 

In  directional  fitting,  therefore,  if  the  scan  ray  intersects  the  virtual  object,  the  intersection  closest 
to  the  instrument  is  chosen  as  the  theoretical  point  P for  the  data  point  P . What  happens  if  the 

scan  ray  of  a data  point  P,  does  not  intersect  the  current  virtual  object?  It  might  be  tempting  to 
reject  such  an  occurrence  as  unrealistic  as  the  point  cloud  was  generated  from  the  real  object.  It 
should  be  kept  in  mind,  however,  that  during  the  fitting  process,  the  virtual  object  will,  in 
general,  not  match  the  actual  object.  Indeed,  establishing  that  match  is  the  purpose  of  the  fitting 
process.  It  is,  therefore,  necessary  to  extend  the  error  definition  to  those  data  points  whose  scan 
rays  miss  the  virtual  object.  The  following  generic  principle  for  a continuous  extension  has  been 

proposed  in  [21],  Here,  the  theoretical  point  P is  chosen  as  a point  on  the  virtual  object  that  is 
closest  to  the  scan  ray  in  terms  of  Euclidean  distance. 


2.2  Sensitivity 

As  we  return  to  the  general  error  function  E (2.1.1),  we  examine  a major  aspect  of  analyzing  the 
results  of  a fitting  procedure.  It  concerns  the  “sensitivities”  of  the  resulting  parameters,  namely, 
their  marginal  rates  of  change  in  response  to  perturbations  of  the  data  coordinates.  Such 
sensitivities  not  only  provide  key  information  about  a fitting  process,  they  also  play  a role  in  the 
estimation  of  variances  and  covariances  of  the  fitted  parameters,  as  will  be  discussed  in  the 
subsequent  section. 
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With  each  set  of  stipulated  data  values  dt,  (pt,  £? , n , the  error  function  E associates  a set 

of  minimizing  parameters.  We  may  thus  consider  these  minimizing  parameters  as  functions  of 
these  data  variables 


(2.2.1)  X{dv(pvO^-,dn,(pn,On) 

Y(dv(p],dl,...,dn,(pn,dn) 

Z(dx,(px,0„...,dn,(pn,On) 

in  a suitable  neighborhood  of  the  actually  measured  values  ,9]"' , i = 1 n.  By 

definition. 

X(0) 

y(0) 

•7(0) 


II 

o 

o 

9;0), 

r\ 

© 

o 

5b 

o 

= Y(di°\<pl0)j 

9,(0), . 

..,dlO),p(O),0ZO)) 

© 

o 

V 

N 

II 

g(0) 

7I  ’ • 

..,di0\(p(0) ,0{O)) , 

expressing  the  desired  results  of  the  fitting  process. 

In  what  follows,  we  will  assume  that  the  error  function  E satisfies  all  necessary  differentiability 
conditions.  We  are  particularly  interested  in  the  derivatives 


dX  dY  dZ  dX  dY  dZ  dX  dY  dZ 

ddt  ddt  dd:  d(pt  d(pt  dcpt  d0[  d9:  d9i 

because  their  values,  for  dt  = d\{)\(pi  = ,9,  = 9'/",  i = 1 n , 

(2.2.2) 


, i = 


dX 


dX 


no)  _ j( o)  /j(u)\ 

a^l  ' '-A  K a* 

| (0)  dY  ✓ .rm  vm.  dY 

dd~ 


1(0)  • 


dX 


(0) 


(0)  _ 


dX  . »rm  /itnu  dx  i I dx 


d<P, 


(0)  _ 

d0  1 d0 


(d,(O),...,0‘O)X  -^-|(0)  = d 1(0),...,C') 


1(0)  • 


dd 


dZ  | (Q)  _ dZ  (Q)  (0). 


d(pt  1 

dZ  i(0)  _ 

dtp,  1 d<p, 


zj(0)i  dY  |(0) 

d(p,  d0,  1 

dZ 


(di°\...,9i{))\  — r = 


3<0)i  | (0) 

d0 


dY_ 

d0, 

dZ 

d0 


-( d[°\...,9™ ) 
(d[°\...,9^) 


represent  the  respective  sensitivities  of  the  parameters  X'0),  Z<0\  Z<0)  to  perturbations  of  the 
data  values  d,01  ,<p}iU  ,9-0) , i = 1,  Implicit  differentiation  will  now  be  used  to  derive 

expressions  for  the  sensitivities  (2.2.2)  from  the  expression  for  the  error  function  E . Indeed,  the 
gradient  of  E with  respect  to  the  parameters, 
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(2.2.3) 


V 


XYZ 


E = 


d 

E(X,  Y,  Z,  dn<p,,0„ 

dx  111 

,.,d  ,(p  , 6 ) 

— E(X,  Y,  Z, 
dY  111 

,.,d  ,(D  ,0  ) 

-^—E(X,  y,  Z,  dv(pvOv. 
dZ 

..,d  ,(p  ,9  ) 

’ n ’ r n ’ n ' 

vanishes  if  the  parameters  X,Y,Z  have  been  minimized  keeping  the  data  variables 
d, , (p, , 6, , ; = l,...,n  , fixed.  As  we  substitute  for  the  parameters  X,  K,  Z their  corresponding 
functions  (2.2.1)  in  the  above  gradient  components,  we  thus  arrive  at  a set  of  derivative  functions 
which  assume  the  value  zero  for  any  stipulation  the  data  variables  d,,  (p, , 6,,  i = 1 . In  other 

words,  the  following  derivative  expressions. 


dE 

dX 

dE 


(X(^, ...,£„),  y ,0„),  z(dr...,on\ dlt...t0n)  = o 


(2.2.4)  ^(X(r/,,...,^),  w Z(</ ,£„),  dv...,Gn)  = 0 

dY 

dE 

teiXid' Y[li 0"l/A(I' = ° 


vanish  identically.  Then  so  do  the  derivatives  of  these  expressions  with  respect  to  data  variables 
- stated  below  only  for  variable  d,  for  brevity  -- , so  that  by  the  Chain  Rule: 


d2E  dX  | d2E  dY  | d2E  dZ  | d2E  = 

dX2  dd , dXdY  dd,  dXdZ  dd,  dddX 

d2E  dX  d2E  dY  | d2E  dZ  | d2E  _ 

dYdX  dd,  dY2  dd , dYdZ  dd,  dd,dY  ~ 

d2E  dX  d2E  dY  d2E  dZ  d2E 

+ + — + = 0. 

dzdx  dd,  dZdY  dd,  dz~  dd,  dd,dz 

Evaluation  using  the  fitted  parameters  X°,Y°,Z°  and  the  actual  data  d" , $>°,  6® , yields  a 
numerical  mxm  linear  system  of  equations  for  the  sensitivities  (2.2.2)  with  respect  to  the  data 
variable  d, . With  the  notations 


ko) 

dx2 ' 


^4(X(0),y<0),Z,0),r/1(0, 

dx2  1 


d2  E 
dXdY 


(0) 


y(0),  Z(0),  J,(0),...,^0)),  etc. 

dXdY 


this  linear  system  takes  the  form 
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(2.2.5) 


d2E  | 

(0)  dX  |(0) 

d2E  | 

o 

QJ 

<0,  , d2E 

hot  dZ 

1(0)  _ 

d2E 

dx2' 

dd) 

dxdy' 

dd,1 

dxdz 

• dd, 

dd,dX 

d2E 

ho)  dX  i(0) 

d2E  |(o)  dY  |(0)  , d2E  \ 

(0)  dZ  | 

(0)  _ 

d2E  | 

dYdX 

dd,  ' 

dY2' 

dd,  1 

dYdZ  1 

dd} 

dd,dY' 

d2E 

1 (0)  dX  |(o)  , 

d2E 

1(0) 

1(0,  . d'-E  1 

(0)  dZ  | 

(0)  _ 

d2E 

dzdx 

dd, 

dZdY 

1 dd, 

dz2 ' 

dd) 

dd,dZ 

The  matrix  of  this  linear  system  may  be  stated  in  terms  of  the  Hessian 


d2E 

d2E 

d2E 

dx2 

dXdY 

dxdz 

d2E 

d2E 

d2E 

dYdX 

dY2 

drdz 

d~E 

d2E 

d2E 

dzdx 

dZdY 

dz2 

of  the  error  function  E . The  linear  system  (2.2.5)  may  thus  be  written  as 


(2.2.7) 


H E 

11  XYZ^\ 


(0) 


~dx 

(0) 

dd, 

dY 

(0) 

dd, 

dz 

(0) 

dd, 

dd. 


V E 

v xkz-M 


(0) 


where  again  the  symbol  l0)  is  meant  to  indicate  the.  — a posteriori  substitution  by 

X Y Z""  and  the  actual  data  values.  The  resulting  Hessian  matrix  is  positive  definite,  and 
therefore  nonsingular  at  any  locally  unique  minimum  of  the  error  function  E . The  linear  system 
is  then  solvable  and  yields  the  values  of  the  sensitivities  (2.2.2)  with  respect  to  d,  and. 
analogously,  for  the  remaining  data  variables,  based  on  the  same  Hessian  matrix. 


Note  that  implicit  differentiation  can  be  used  to  determine  higher  order  sensitivities  such  as 


d2X  d2X  d2X  d2X  d2X  d2X 

dd~  dd,d(pr  ddid0J  d(p~  d(p,d9  d6 


The  corresponding  linear  systems  are  based  on  the  same  Hessian  matrix  as  in  (2.2.7)  but  use 
different  right  hand  sides. 
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2.3  Noise  Propagation 


In,  general,  data  variables  d,,(p,,6,  and  the  parameters  X,Y,  Z will  be  considered  random 
variables  with  expected  values  d^  and  X(0),  T(0),  Z(0),  respectively.  In  many 

applications,  however,  some  data  variables  of  the  error  function  will  be  “control  variables”  or 
“design  variables”  which  are  given  and  thus  not  random.  When  fitting  scanned  objects,  it  is 
commonly  assumed  that  the  noise  in  range  measurements  di  dominates  the  noise  in  bearings. 
Consequently,  only  the  range  coordinates  d,  are  considered  random,  while  the  bearing  angles 
(p,  and  0,  are  specified  control  variables.  For  scanning  instruments,  it  is  generally  safe  to 
assume  that  range  variables  dt  are  independent  of  each  other.  The  following  exposition  will  be 
based  on  these  assumptions. 

The  sensitivities  described  in  the  previous  section  will  be  instrumental  in  assessing  the  effects  of 
data  noise  on  fitted  parameters.  The  well  known  “Error  Propagation  Formula”  provides  first 
order  estimates  of  the  variances  (see  GUM  [23],  Chapter  5) 


(2.3.1 


var(X)  = ^ 

/=1 

n 

var(T)  = ^ 

i= 1 

n 

var(Z)  = ^ 


i=i 


dX  | 
dd) 

dY_\ 
dd,  1 

dz 
dd . 


(0) 


(0) 


(0) 


var (di ) 
var  (d,) 

) 

var(c/,) 


Similarly,  one  has  for  the  covariances  [24], 


i=i 


(2.3.2) 


cov(T,Z)  = y 

/'  = 1 

n 

CO  v(Z,X)  = y 


/=1 


r dx 

o 

I(0) 

{dd, 

1 dd, 

1 

(0)  dz 

(0) 

dd, 

dd) 

— i 

Q J 

N 

U0)dX_ 

(0) 

dd, 

var  (d,)  = cov(T,  X ) 


var {d, ) = cov(X,Z) 


Note,  that  even  as  the  measured  quantities  d.0)  are  independent,  the  fitted  parameters 
X(0),  T(0),  Z101  will  still  be  correlated. 

In  most  applications,  the  condition  of  homoscedacity  is  supposed  to  hold:  Tat  is,  the  variances 
have  the  same  value  within  a class  of  measurements,  such  as  the  class  d of  range  measurements, 
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var(d)  = var(c/ ),  i = 1 ,...n  . 


For  the  special  cases  of  LS  and  NLS  regression,  the  individual  errors  in  our  scenario  would  take 
the  special  form 


The  matrix  of  variances  and  covariances  can  then  be  approximated  in  an  elegantly  simple  fashion 
under  homoscedacity  as  set  forth  in  the  general  literature,  e.g.  [5-7,  25,  26].  Unfortunately,  the 
error  functions  considered  here,  in  particular,  the  orthogonal  error  function  (Chapter  4)  and  a 
portion  of  the  directional  error  function  (Chapter  3)  do  not  fall  into  this  regression  category 
because  the  required  separation  of  the  random  variables  from  the  control  portion  cannot  be 
achieved.  It  is  for  this  reason,  that  the  more  general  approach  described  in  Sections  2.2  and  2.3 
had  to  be  adopted. 


3 Directional  Fitting  of  Spheres 
3.1  Directional  Errors 

Introducing  the  trigonometric  quantities  ^,77,,^,  the  Cartesian  coordinates  xi,yi,zi  of  data 
points  will  be  expressed  in  the  form 

(3.1.1)  x,  = di  cos  (pt  cos  6i  = = d;  sin  <p.  cos  0t  = dirji , z,  = d,  sin  <9  = diqi , 

where  + T]~  + q~  = 1 . The  vector  (^  j]i , ) represents  the  direction  of  the  scan  ray  along 
which  the  data  point  P was  acquired.  Next  we  introduce  the  quantities: 

p^X^+Y^+Z^, 

(3.1.2)  qt2  = X2  + Y2  + Z2  - pi*  , 

7 _ ->  ? 

s-  = R--q~. 

Figures  1 and  2 illustrate  the  geometric  meaning  of  these  quantities. 

If  the  scan  ray  of  data  point  P intersects  the  virtual  sphere  centered  at  C = [X  f Z],  we 
associate  with  P the  midpoint  P of  the  resulting  chord.  The  quantity  pt  =|  P' ' then 
represents  the  distance  of  that  “mid-chord"  point  from  the  instrument  at  the  origin  O =[000], 
Similarly,  q ( = \\P'  -C ||  > 0 represents  the  distance  of  the  mid-chord  from  the  sphere  center.  It 
follows  that 

q,<R 
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is  the  condition  for  “true”,  that  is,  non-tangential  intersection.  The  quantities  p:  and  qi  are  side 
lengths  of  the  right  triangle  AOP'C , with  its  right  angle  at  P' . Pythagoras  thus  yields  the 
relation  (3.1.2)  between  p~  and  qi  . The  triangle  AP  P C , where  the  theoretical  point  P 
marks  the  entry  point  into  the  sphere,  also  has  a right  angle  at  P' . The  quantity  y thus 
represents  the  length  of  the  half-chord  that  needs  to  be  subtracted  from  the  distance  p ( = ||P'||  to 


arrive  at  the  desired  distance 


of  the  theoretical  point  from  the  origin. 


As  a result,  the 


directional  error  of  the  data  point  P is  given  by 


(3.1.3)  f]  - Pi  ~ y - di  (-  “interior  error”  of  P if  qi  < R) . 


Figure  1.  Geometrical  interpretation  of  the  directional  error  function  if  a scan  ray  intersects 
the  sphere  surface.  P,  (marked  by  dark  dot)  is  the  experimental  point,  light  dots  mark  the 

theoretical  point  on  the  sphere  surface  P and  the  mid-chord  point  P . The  length  of  the 

bold  line  segment  measures  the  error  f]  defined  by  (3.1.3). 


On  the  other  hand,  if 


<7,  ^ R , 

then  the  scan  ray  of  the  data  point  P fails  to  truly  intersect  the  virtual  sphere.  In  that  case,  we 
follow  the  general  extension  principle  set  forth  in  Section  2.1  aqnd  determine  the  theoretical 
point  as  that  point  on  the  virtual  sphere  which  is  closest  to  the  scan  ray.  The  line  segment  which 
represents  the  shortest  distance  between  the  sphere  and  the  scan  ray  has  to  be  orthogonal  to  both 
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the  sphere  and  the  scan  ray.  The  line  segment  thus  has  to  be  part  of  a line  through  the  center  of 
the  sphere,  and  also  meet  the  scan  ray  at  a point  P'  at  a right  angle.  This  defines  again  the  right 

triangle  AOP'C , which  we  encountered  before,  and  whose  side  lengths  are  again  pt  and  qt 

(Fig. 2).  The  desired  theoretical  point  P thus  lies  on  the  side  [ P'O  ] at  distance  R from  center 

and  at  distance  q,  - R from  f*' . The  triangle  AP  P C is  also  a right  triangle,  has  side  lengths 


\P  P II  = p - d and 

/ /II  r i i 


P'-P 


= qt-  R . The  length  of  its  hypotenuse 


P -P 


thus 


represents  the  error  of  the  data  point  P : 


(3.1.4) 


8i 


yjiP,  ~ df  + {q,  - Rf 


(=  "exterior  error”  of  Pt  if  qt  > R ). 


Figure  2.  Geometrical  interpretation  of  the  directional  error  function  if  a scan  ray  does  not 
intersect  the  sphere  surface.  Pi  (marked  by  dark  dot)  is  the  experimental  point,  light  dots 

mark  the  theoretical  point  on  the  sphere  surface  P and  the  point  P on  a scan  ray  which  is 
the  closest  one  to  the  sphere  center  C.  The  length  of  the  bold  line  segment  measures  the 
error  g_  defined  by  (3.1.4). 

If  r/(  - R = 0 then  ft  — , so  that  the  combined  error  function  will  be  continuous.  While  the 

error  expression  g,  is  everywhere  twice  continuously  differentiable,  the  error  expression  fj  fails 
to  be  so  if  and  only  if  st : = 0 , — the  case  of  tangential  intersection  — , where  its  gradient  with 
respect  to  the  parameters  X,  Y.  Z is  infinite.  In  these  cases,  the  resulting  full  error  function 
will  also  not  be  differentiable.  However,  those  points  will  only  amount  to  a closed  set  of 


measure  zero  in  parameter  space.  As  a consequence,  a gradient  based  numerical  minimization 
method  such  as  the  often  relied  upon  “BFGS”  method  [27]  may  still  be  used  [21].  Similarly,  the 
probability  of  the  error  function  E not  being  differentiable  for  the  fitted  parameters 
X<0),  T(0),  Z""  will  be  theoretically  zero. 

We  find  it  convenient,  to  categorize  only  a true  intersections  as  a “hit”.  A tangential  intersection 
is  thus  considered  a “miss”,  along  with  all  cases  in  which  the  virtual  sphere  is  not  met  at  all. 
Accordingly,  we  divide  the  indices  / into  two  sets: 

hit  - [i : qt  < /?)  and  Ext  - [j  ■ q ; > /?}  . 

The  combined  error  function  then  takes  the  form 


(3.1.5) 


EJn,  = Edrc, (X,  y,  z,  d, , n , ....  d, , 0,  ,0„)  = Y ft  + £ 


* je  Ext  ^ i 


3.2  Derivatives  for  the  Directional  Error  Function 

In  this  section,  we  list  formulas  for  the  gradients  and  the  Hessians  of  the  individual  error 
functions  f{  and  gl  with  respect  to  the  parameters  X,  T,  Z,  along  with  the  second  derivatives 
with  respect  to  both  these  parameters  and  the  data  variables  dt . , (pt . , 0t . Gradients  support 

optimization  methods  and  are  the  first  step  towards  determining  the  above  second  derivatives, 
which  are  needed  for  the  computation  of  the  sensitivities  and  variances  described  in  Sections  2.3- 
4.  Derivation  of  these  formulas  is  provided  in  the  Appendix  as  referenced. 

In  terms  of  the  individual  errors  /,  and  gt  , the  gradients  and  Hessians  of  the  directional  error 
function  are: 


^ XYZ^drct  L/.U  + ZhieExy  XYZS, 

H Xr/Ejrc,  = Z 'jieln^XYzfi  + Z,e  ExtH^gi  ' 

Gradients  are  considered  column  vectors  and  are,  in  general,  linear  combinations  of  the  vectors 


"x" 

A 

(3.2.1) 

U = 

Y 

and  At  - 

n, 

Z 

Jdi  _ 

For  instance  (see  (A. 2. 7-8)  in  Appendix  A), 
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(3.2.2)  -Vxrzf>  = 


Ilu 

m+dj  i 

i ’ 

iv  p 2 = 

~ v XYZ  6 i 

r._*i 

u + 

[*>-<) 

s i 

i 

2 

l <h) 

l Pi  J 

Similarly,  the  Hessian  matrices  are  linear  combinations  of  the  following  four  symmetric 
matrices: 


1 0 0 


(3.2.3) 


/ = 


0 


(Identity), 


0 0 1 


~x~ 

'x2 

XY 

xz~ 

UUT  = 

Y 

\X 

II 

"no1 

X, 

YX 

Y2 

YZ 

z 

zx 

ZY 

Z2 

~X~ 

'Xl+lX  Xr,,+£Y  X<r,+£Z‘ 

UAj  +AiUt  = 

Y 

[£  p Cih 

Pi 

X 

x 

N 

II 

Y£,  + r 1,X  Yrj,  + n,Y  Yg+ri.Z 

z 

_Gi  _ 

_Z^+g,X  Zrj:  + g.Y  Zg,+f,Z_ 

i 

I 

1 

P& 

£ r 

44r  = 

Pi 

II 

Cr» 

Cr 

p& 

1 

p: 

P&i 

_bi  _ 

sa 

2 

Thus,  (see  (A. 2.1 1-12)) 


(3.2.4) 


-HXYZf*  = 


+ 


P,  ~ d, 


\ 

/ + 

Pi  ~d, 

UUT- 

\p,(p-d,)  l) 

) 

l s>  ) 

l si  si  J 

(UAT  + a,ut  ; 


(pf  -sfXPj-dj)  2(pi-si) 


A.  A 
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f 


Hxyz  gf  = 


R 


(3.2.5) 


1- 

V ch  J 


I +-^UUT  -^(UAj  +A,Ut) 

(i  q\ 


+R\ 


x2  + y2  + z2^ 


q, 


A A 


Evaluated  for  optimal  parameters  X(0),  Y{0\  Z(0)  and  actual  data  points  Pl  = [d]"'  (p\Q]  6>" ' ] , 
these  Hessian  matrices  support  the  left  hand  side  of  linear  systems  (2.2.7)  for  the  corresponding 
sensitivities.  For  the  right-hand  sides  of  those  systems,  we  have  for  the  range  variables  di  , 


v — e = y v ^1  + y v ^ 

dd  drec  ^ iel,u  XYZ  dd 


dd. 


where  (see  (A. 3. 1 -2)) 


(3.2.6) 


-V 


a/,: 


2 xr/  dd, 


— -U  - Pl — — A 


-V  XYZ  — = ~A 

2 dd, 


For  the  bearing  variables  (p  , Q , we  obtain  the  individual  derivatives  in 


vx„—Ed  = y vm'dl~  + Y vxrz&- 

XYZ  drec  Jielnt  XYZ  ^ ^ / w /e Ex/  AY/ 


9ft 


9<P, 


v —e  = y v ^ 

v XYZ  00  Z-f/etor  V *VZ 


9ft 


2 + y v ^ 


do. 


by  multiplying  the  vectors 


-n, 

(3.2.7) 

,*4  = 

d(P, 

p 

. a4  = 

9? 

Pi 

0 

7, 

where  a,  = -cos# .sin 0t /?,  = -sin #.  sin <9, , y;.  = cos  <9, , by  matrices  /7(/.2)  and  /7(g,2)  as 
follows: 


-Hi 

(3.2.8) 

’v  ^ =r(f2) 

v xrz  -x  1 ’ 

2 9# 

i 

yy  o 

i 

1 y 9S,2 
’ 2 XYZdcpi 

= eu,2) 

p 

0 
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(3.2.9) 


2 


XYZ 


<*i 

a, 

r,(f,2) 

*v  = r (t>2) 

2 XYZ  d6,  iK8i 

Pi 

Yi_ 

Yi_ 

Here  the  two  pre-multiplying  matrices  (see  (A. 3. 8), (A. 3. 1 1))  are  given  by 


(3.2.10) 


-r,(fr)  = 


Pt(P,-d,) 


v5- 


+ (Pi -Si) 


f ,\ 

Pi  - d, 


UUT+(p,-s,) 

/, 


'(Pi+SiXpi-dj)  2a 


AUt 


(3.2.11) 


r,(g;)  = 


f D > 

RPi 

3 

UUT  + R 

l Pi  J 

V 

f x2  + y2  + z: 

3 

Pi 


\ 

a,ut  + 

( n \ 

Rp'-d 

) 

l Pi  J 

Because  of  the  following  orthogonality  relations. 


(3.2.12) 

4 

£ 

= 0 and  A] 

P 

0 

_r,_ 

the  above  matrices  may  be  replaced  in  (3.2.8)  and  (3.2.9),  respectively,  by  the  following 
symmetrized  versions: 


-r5(/2)  = 


'l  Piip-df 


(3.2.13) 


+ (Pi~Si) 


P,  ~ d ' 


UUT  + (p,  - 

/, 


(Pi+SjXp.-dj)  2 
S 


(AtUr  + UAj ) 


’<  J 


(3.2.14)  -Hig;) 


_RPi_ 

V Pi  ) 


UUT  + R 


( X2  + Y2  +Z2^ 


(AUr  + UA’)  + 


RPi 


d. 


\ p, 


i . 
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4 Orthogonal  Fitting  of  Spheres 


Here,  the  theoretical  point  PI  , that  is,  the  point  on  the  sphere  which  is  closest  to  the  data  point 
P = [x;  yi  z(],  defines  the  individual  error 


P-P 


= \\C-PI\\-R  = w.  - R 


with  respect  to  the  sphere  center  C and  the  radius  R . We  thus  represent  the  orthogonal  fitting 
approach  by  the  full  error  function 


(4.1.1)  = 2>,2  = 2>,  - R)2  = K ~ 2 Rw,  + R2  . 

(=1  (=1 

Consistent  with  the  generation  of  point  clouds  by  scanning  from  a single  instrument  location,  and 
as  discussed  before,  the  underlying  coordinate  frames  are  again  considered  polar  with  the 
instrument  location  at  the  origin: 


xi  = di  cos  (pt  cos  O'  = d'^ , yi  = di  sin  <pt  cos  (9,  = djl , z,  = dt  sin  0t  - d£t . 

For  an  analytic  discussion  of  the  orthogonal  error  function  in  terms  of  Cartesian  data  see  [11]. 
With  the  definition  (3. 1 .2)  of  the  auxiliary  quantitity  pt , we  have 

(4.1.2)  wf  = (X-xi)2+(Y-yi)2+(Z-zi)2  = (X2 +Y2  + Z2)  - 2d,p,  + d2  . 

A key  vector,  in  which  gradients  and  Hessians  of  the  individual  orthogonal  error  squares  hf  may 
be  expressed,  is  given  by 


(4.1.3) 


£ 

zn 

w,= 

Y-yt 

= 

Y 

~d, 

n, 

_z-z,_ 

Z 

= U - d,A, 


as  W/W  = w2  . Also  (see  Appendix  (A. 4. 3)  and  (A. 4. 5)), 


(4.1.4) 


-VxyzK 


R 


v 


W 


and 
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(4.1.5) 


\Hxrzhf  = 


R 


V vv<  J 


R T 
I + —WWr 

3 i i 

w. 


For  the  derivatives  which  define  the  right-hand  sides  of  the  linear  system  (2.2.7),  we  have  first: 


XYZ 


dd, 


e = y v 

t-'orth  V 


XYZ 


1=1 


av 

dd, 


with  (see  (A. 5. 2)) 
(4.1.6) 


ln  dhf  R/ 

~ ^ XYZ  — 3 ( Pi  dj  )Wi 

2 ad;  w; 


v w,v 


A . 


Again,  the  corresponding  mixed  derivatives  with  respect  to  the  bearing  variables  (pt , 6 

are  multiples  of  the  vectors  -^—A  and  —A  , defined  in  (3.2.6).  Their  common  multiplier  is  the 

d(pt  ' dO, 

matrix  (see  (A, 5. 4)) 


-r,(tf)  = -4 


R 


V wij 


I + Ml wut 

w 


(4.1.7) 

which,  in  analogy  to  (3.2.1 1),  may  be  replaced  also  by  its  symmetrized  form, 

(4.1.8) 


~r'(h;)=~d 


d-«d 

wu 


/ + ELwyvJ  . 

w; 


Thus 


-Oi 

-o, 

(4.1.9) 

]vXYzdh‘  = n/*2) 

2 dtp, 

£ 

= Mtf) 

£ 

0 

0 

and 


(4.1.10) 


-vXYZ  — = T(h2) 
2 dO 


V 

V 

Pi 

= rtf) 

Pi 

Yi_ 

Yi_ 
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Appendix  A:  Determination  of  Derivative  Formulas  Used  for  Calculating  Sensitivities 

Here,  we  provide  step-by-step  developments  of  the  derivative  formulas  referred  to  in  Chapters  3 
and  4 for  the  purpose  of  determining  the  parameter  sensitivities  for  directional  and  orthogonal 
fitting.  In  Section  A. 2,  the  gradient  Vxyz£/rec  and  Hessian  H XYZEdrec  of  the  directional  error 

function  are  at  issue.  The  Hessian  provides  the  matrix  for  the  corresponding  linear  system 
(2.2.7).  Also  for  the  directional  error  function,  the  derivatives  of  both  parameters  and  data 
variables,  are  derived  in  Section  A. 3,  furnishing  the  right  hand  sides  of  these  systems.  Finally, 
Section  A. 4 provides  the  analogous  information  in  the  case  of  orthogonal  fitting. 


A.l  General  Formulas 


In  what  follows,  the  calculation  of  gradients  Vxyz  and  Hessians  H XYZ  will  often  be  based  on  the 
following  straightforward  reformulations  of  product  and  chain  rules: 

(A. 1.1)  V XYZb(a)  = b'(a)V XYZa  , 


aV  XYZb  + bSZ XYZa  , 

^ xyzci  ~ 2aV  XYza  , 


and 


(A.  1.2)  H XYZb(a)  = b"(a)V XYZaVTXYZa  + b\a)H XYZa  , 

H XYZab  = V xYza^TxYz^  "b  V xYz^J\rza  aH xvzb  4*  bH XYZ a , 

H XYZa  ~ XYZa^ X\7,a  + 2 ClH xYZa  ’ 

H xvza  ~ H XYZ xl ci  = ——  V XYZa  4-—  H XYZ a . 

4 n La 


These  formulas  are  straightforward  reformulations  of  Product  and  Chain  Rules,  and  will  not 
always  be  referenced  in  what  follows.  Gradients  of  functions  are  considered  column  vectors. 
Gradients  of  row  vectors  will  result  in  matrices  as  the  elements  ol  a row  vector  will  become 
columns  of  gradients.  For  vectors  U and  At  (3.2.1),  we  thus  find 
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1 

0 

0 

0 

0 

0 

(A. 1.3) 

V UT  - IV  XV  TV  Z 1 = 

y XYZ^  |vXKZyv  y XYZ1  y XYZ^  \ 

0 

1 

0 

= I , 

V,E  4 = 

0 

0 

0 

0 

0 

1 

0 

0 

0 

Formulas  (3. 2. 8-9)  and  (4.1.9-10)  are  based  on  matrices  which  pre-multiply  derivatives  of  At . 
These  matrices  are  the  result  of  applying  the  differential  operator 

d2 

dXdg, 

a2 

MS 

a2 

dZdg 

to  error  functions. 


(A. 1.4) 


r = V VT 

1 , y xrzy  in,?, 


a2 


ara^ 

a2 


a2 


dxd£  dXdr/, 

a2  a2 

dYdrj, 

a2 


aza£  aza/7, 


A.2  Gradients  and  Hessians  of  the  Directional  Error  Functions 

Recall  the  directional  error  function 


E - = V f2  + V p 

Jrc7  Z— lieU  J 1 —it  l'  * ' 


with  individual  errors  (3. 1.3-4), 


fi  = Pi-  s,  ~ d,  ’ <?,  = -df  + ^-RY 


based  on  the  auxiliary  quantities  p, , qt , s,  (3. 1 .2)  and  the  direction  cosines  (3.1.1) 


<g(  = cos  (pt  cos  0: , rji  = sin  (pt  cos  6? , ^ = sin  Gt 


Using  (A.  1.1)  where  indicated,  we  note: 


(A.2.1)  = Vxlz(Xg  +Z77,.  +Zffi)  = 4 , V XYZ p2  = 2p,4  , 


V*^,2  = V,^(X2  + T2  + Z2)  - VXY7p2  = 2[U  - pA\  , 

V*^,  = ^ xyz^  = \~S/  xyz^'  =—  - P,  4 1 , 

2 <1,  q, 
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(A.2.2)  VXyzsi2  = Vxyz(r:  ~qf)  = ~v xncll  =~2[u  - P,A]  , 

V XYZS,  = ^XYZyR  =\— VxYZSi2=-—  lU  ~ PiAi\  * 

2 Jf  5, 

= T ^ xYzsi  =Mu  - PiA]  ■ 

s,  s,  s' 


All  Hessians  H XYZ  determined  in  this  section  will  be  linear  combinations  of  the  four  matrices 
(3.2.3).  Again  we  begin  with  the  auxiliary  quantities: 


(A. 2.3)  HXYZpi  = 0 , Hxyzp;  = 2W  XYZpiVrXYZpi  = 2A,A^  by  (A.  1.2), 

Hx„  q;  = Hxrz(X2  + y 2 + Z2)  - Hxrzp2  = 2(/  -A,AJ ) , 

Hxrzsf  = - Hxrzq, 2 = -2 (/  -44r)  , 


From  Hessians  of  q] , s,  we  pass  to  Hessians  of  qt  , si , using  (A.  1.2) 


Thus 


-1  1 


H XYZCli  ~ ^ XYZ  V qi  A 2 ^ XYztfi  ^XYZCh  + ~ ^XYZ^i 

4 q.  2 q, 

= -\(U-  p,4  )(U  - P,A,  f + -(/-  A, A2 ) 

qi  <k 

~ uuT  + 4 (t/4  + AU ) - 4 44r  + — (/  - 44r )■ 

q,  q\  q\  q> 


(A. 2.4)  HXYZq,  =-I-\uUT  + ^f(UAi+AiU)- 

q,  q\  q\ 


C 2 

P, 


+ 


q\  q 


A A 


i y 


Concerning  the  last  term,  note 


2 t 2 2 2 

Pi  i Pi  + _ x ' + ^ ^ 

2 + 2 3 

q,  q,  q,  q, 
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Similarly,  by  (A.  1 .2), 


ww  tw  / 2 ^ *■  v-?  2\~iT  2 ^ 1 * > 

™ XYZSi  ~ **xYZ\Si  ~ ~~A  XYZSi  XYZS  i ~Z  H XYZS i 

4 s;  2 s, 

= __!_([/_  pA)(V  _ pAj  - — (/  - 44r ) 

5.  .S'. 


1 


= -- j UUT  + (UA,  + A,U)--^-A,Aj  --(I  - 4,4 

S-  S';  S';  S. 


Pi 


Thus 


(A.2.5) 


H XYZS I 


— + +4-t/)  + 

S;  S'.  S'; 


2 \ 


\Si  Si  J 


4;  A7 


Concerning  the  last  term,  note 


p;  _ s;  - p;  _ R~  -q;  - p-  _ R~  ~(X~  + Y +Z~)  _ X +Y~  + Z - R~ 


S;  s] 


Vi 


<h 


<h 


<h 


From  the  above,  we  derive  derivative  expressions  involving  the  errors  f,  and  g, : 


(A. 2.6) 


Vxrzf,  = VxyzP,  - vxrzsi  = A + —{U  - pA, .) 

*,■ 

= lU.£i^lLA=lU.A±^Ai 

S S;  S.  S. 


(A.2.7) 


Atvzf.2  = fSxrzf,  = 4u  - f-,J  +l/  '4 
2 5 5 


For  the  external  portion  of  the  error  function,  we  find 

V xyz 8 i~  = 2(P,  ~dPVXYzP,  +2(^,  -R)V xyzP,  = 2( Pi  - d,  )U  + —q,-—{u  -p,A,) 

Qi 

2(p i — R)  r r , 2(p-d,)q,-2(q-R)p, 

q,  q, 

or 

(A.2.8)  A*>?S'2  ={{--)U  + {—~d)  4 

2 q,  q, 
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and  by  (A.  1.1), 


V e = V ,/p  2 = -—V  g 2 — — 

v XYZ  6 / y XYZ\  6 , - y XYZSi 

2 g,  g, 


(A.2.9) 


l <h  <h  j 


d,g, 


■{(q,  -R)U  + (R p,  ~diqi  )Ai)  = — — 11 U + — ^ — ^^-A: 


- q‘~Ry  | RP<-di(k 
dig  i d,g, 


Moving  to  the  second  derivatives,  we  find 


(A.2.10)  Hxrzf,=-Hxyzs,  =-I + \(JUT  ~^(UA[ + A,Ut)  + 

S:  s ; 5; 


Pi 


2 \ 


m: 


Next,  we  introduce  the  matrix 


<W,  = V,K/VrXK/,  = —{V  -{p-sM)(U -(P,-S,)A,J 

Si 

= \uUT  -^-^(t/47'  +A1UT)+(P‘~1S‘)  a, a\ 

s'  s'  s' 


Also  by  (A2.10), 


fi  ' • firrrrT  fi  Pi  / j i aT  , ajjTs  , f,(Pi  sf)  T 


f,HXYZfi  - —I  + ^jUU  - 2-Ly- (JJA‘  + AjU  ) + 


A A 


3 r 1 


By  (A.  1.2),  Hx„fr  = 2Gmf,  +2f,Hxrzf,  . Thus 


\nmfi=Z-I  + 

2 s 


+ 


fi 


- ^ 


UU‘ 


A 9 


Pi  ~ si  , fiPi 


v 


+ 


2 s \ 


(Pi-Sj)2  ! fj(pf-sf) 


AA  . 


9 


(t/47-  + 4^) 


Expressing  f]  = pi  - y -di  yields 


+ 


/,  _ s,  + /,  _ p,  ~ dt 


1 2 

S 5- 

/ / / 


Pi  - s<  | ./>,  = ip,  -s,)Sj  +./>,  = /Vy,  - ^ + /V  - PA,  - M = Pfpt~d,)  _1 
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(Pi-Sj)2  | fiipf-sf)  = pfSj  - 2 Pis;  +S*  + p-  - p,s;  - siP;  + ,v;  - p;d,  + d/; 

3 ' 3 3 

5*;  5;  S; 


_ ( P',  - PA1, ) - ( p/;  - d,X ) - 2P ,-V  + 2T 


( p]  - s] )( Pi  -di)-2(pi-  s,.  )sf 


From  these  relations,  we  find  the  alternate  expression 

x f 


1 


nXYZf,2  = 


(A. 2.1 1) 


Pi  ~d, 


v */ 
( 


1 + 


J 


Pi~d, 


UU 


V ‘ J 


Pi(pj-d,)  1 

V 5 1 s 


(UAj  +A,UT) 


i y 


+ 


f AX  -s;)(p,~d,)  | 2(pi-si) 


A A 


which  establishes  (3.2.4). 

We  move  to  the  external  part  gt  of  the  error  function.  The  Hessians  H XYZpt  di'  H XYZd  i’  H XYZ^ 
vanish.  Thus 

Hxyzg;  = Hxyz((  p.  - dt  y + (q, - Rf ) = H XYZ  ( p]  - 2 pidi  + d2  + q2  -2  qtR  + R: 

= HXYZ(pf+qf-2Rqi)  = H XYZ(X 2 +Y2  + Z2  - 2 Rq.)  = 21-  2 RH XYZq, 


and  by  (A. 2. 4), 


(A.2.12) 


Hmgt  = 


R 

V dj 


R 


I +d\UUT  ~-^(UAt  +AUt)  + R 


f X2  + Y2  + Z2  ^ 


<fi 


AA 


This  establishes  (3.2.5). 


A.3  Mixed  Derivatives  of  the  Directional  Error  Function 

The  right-and-sides  of  the  system  of  linear  equations  (2.2.7)  are  at  issue.  They  require  the 
negatives  of  mixed  derivatives  of  the  form 


— V E 

9*  yXYZ^rc, 
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where  * indicates  a data  variable  of  the  error  function.  We  first  consider  the  data  variable  di . 
Note,  in  this  context. 


V 


XYZ 


dd, 


= Vvv7(-l)  = 0,  v 


XYZ 


dd 


= V 


XYZ 


V 


^ XYz(  2 f /)  ~ 2^S XYZ  Pi  ^ XYZSi) 


Thus  by  (A.2.1), 


(A.3.1) 


-V. 


b XYZdd 


-A,  -—(17  - p A ) = — -U  + (— - 1)4) 

S,  S;  S 


(compare  (3.2.6).  Similarly, 


^ xyz  0^  ~ ^ xyz  {(Pi  +{q,  R)  )—  2(Vxkz(/?/  ) — ~^xyzP, 


so  that  (compare  (3.2.6)) 


(A.3.2) 


lv  M 

2 *rz  dd. 


= -/l  . 


As  pointed  out  in  Section  3.2,  the  calculation  of  the  corresponding  derivatives  with  respect  to  (pt 
and  6t  will  be  based  on  the  matrices  generated  by  the  differential  operator  (A.  1.4) 

r = . . 

In  order  to  apply  this  operator  to  the  individual  errors  f2 , g2  , we  first  apply  it  to  the  auxiliary 
quantities  (3. 1 .2).  For  a first  stage,  we  note: 

(A.3.3)  V^Pi . = VllJX£+Yrli+Zgi)  = [x  Y z]  = UT , 


= 2py^,p,  = 2 PiU‘ 


%n^?=Vlrt,JX2+Y2+Z2  - pf)  = -V'^/v  = ~2PiU 


7 2 


a = .[q2  = - * q2  — ——UT 

di 


A V4, 


V*  — = *-V7  a = —UT 

v £-7;C,  2 v £ -She -Hi  3 ^ 

d,  d,  d, 
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(A. 3.4) 


= !--q;)  = ^iq;=2PV 

S/T-  5 = V e a/7  = -^Lv^  s2=—Ut 
Z.n.s,  ' v < ->  FT  l,n.s,  ' „ 

z v-v;  -v- 
Vr  -1  = Lyr  ? - -^-UT 

S:  S2  «**'  " 4 ’ 


and  consequently. 


(A.3.5)  Vlr]i?ipi-si)  = 


( \ 
Pi 


U 


r=-h^uT 


v A y 


v k?,  ( Pi  - s> y 2 = 2( p,  - ) v L,?, ( p>  -*.-)= _2  ( /?'  v ~v' } ^ 7 


For  the  next  stage,  in  view  of  the  above,  (A.2.1-2),  and  the  Product  Rule  yield 


r ( Pi -s.)  = V xyz V ^ ( P,  ~ s, ) = V 


xrz 


-P^UT 
v y 

y 


(VxYztP.-s,)^'  - (Pi-Si) 


A - 


-(u-pA) 


v s. 


V — 

v xrz 

v *.-y 
^ y 

UT-(p,-s,) 


Ut-^VxyzUt 


JJ 


\s. 


T (U-Pi4) 


UT-^1 


u + 


y \ 

> 

( 

1-^L 

4 

uT- 

V 

p1z!lvJp,-si)p,/. 


Ur-J^l 


iU  + h^LA 


V 


.2  I 


t/r  + 


Pi  - V.  U ! ( P<~S,)P,  A 


uT  - Pah.  i 


so  that  collecting  by  matrices. 


ri(P,-s,)  = 
yields 
(A.3.6) 


1 P,  - *,• 


v si  si  y 


/7t/r  + 


f Pi~  si  , (P,-S,)P  ^ 


\ 


AVT_PiZlL, 


y 


FiiPi-Si)  = 


y _ 3 y 

_ Pj_ 

UUT  + 

\ s<  ) 

V 

(P  +si)(Pj-si) 


AUt  - — — — / 
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Next,  by  (A. 3. 5), 


r,  ( Pi  - Si  )2  = V xkVJ(?7i?i  ( p,  - x,  )2  = V 


xrz 


,2  X 


2— — 


x. 


£/r  = -2V 


XKZ 


= -2-(Vxk(A -^)2^7  -2(pi-siy 


f 


(p,~s,y  vt 


V 


V — 

v XYZ 

V ■Sy 


UT  —VxyzUt. 


Referring  back  to  (A.2.1-2),  we  find 


^ XKZ  ( Pi  Si  ^ ^ XYZ  Pi  ^ XYZ  Si  ^ 


V 


P(U  - p,A,)\  = -u --£—±4 

j 5,  x,. 


v ( Pi  - si )2  = 2( Pi  - x, ) V ( Pi  - Si)  = 2( /?,  - x, ) 


2 (Pi-sjy  2{pi-sif  A 


V' 


x, 


; y 


Substituting  the  latter  expression  in  the  above  expression  for  /",(/?,  -x,)2,  as  well  as  by  (A. 2. 2) 
and  (A.  1 .3),  we  arrive  at 


X-r,(P,-sf  =-! 

2 5; 


2 X 


2<P,  -■S)^  2(Pi-Sj)  A 


V 


- (Pi~  s,y 


>-M) 

vj; 


(P.-Sj)-  j 


y 


x, 


= + 


+ 


X-  X-  J 

(A-^)2.i/  + (/?,-^)2A  ^ W (P^yY 


-A.  U‘ 


/, 


so  that  collecting  by  matrices  yields 


-r.ip.-yy  = + 


2(  Pi  ~ x, ) (/?,-x() 


2 \ 


t/l/3 


+ 


UPi-s^  Pi(Pi-Si) 
; 1 ; 


2 \ 


V */ 


AUT  — ^-/ 


’/  y 


X. 
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and 


(A.3.7) 


-r,(P,-s,)2  = 


6 2 . 2 
- P,  + A 

2 

/ 

t/£/r  + 

l A J 

V 

(P,  -Si)  (Pi  + 2s i) 


AnT  ~ih — £il_  / 


As  f2  ={pi-slf  +2di{p,~sl)  + d; , we  have  — P,(f,~)  = s, ):  -ttr^p,-  si ) 

so  that  by  (A. 3. 6-7), 


/ 2,2  y \\ 


pi  + a 


d. 


\ 


s. 


Pi_ 

V J) 


uu1  + 


(p.-sfiP'+ls,)  d (Pi  + si)(p,  -s,) 


A.U‘ 


+ 


f (Pi -Sj)2  d (pi-si)^ 


\ 


/. 


Now, 


2 2 
pi  + a~ 

9 
iS* ; 


~d. 


f ^ 

Pi 


;r  + 5;  + </,;>,  1 p,  (;>,-<) 

.S'2  .S'.  .S',3 


(P/  ~N):(A  +2N)  _ d (Pi  + A'XP,  -5,.)  = ( A - A )2  ( /?,■  +*,-)  + ( P,  - A A A - d,  ( pt  - 5,.  )(  P,  + ) 


( _s  )(P,  + A )( Pi  ~ s, ) + ( P,  ~ si  )s,  ~ dj  ( Pl  + -S', ) 

53 


= (P,  ~ A) 


(p,  + a)(p,  - a ~di)  + (Pi  - a)a 


= ( _ _y  } ( Pi  + A )( P,  ~ O ~ ( Pi  + A A + ( P,  - A ) A 

s.3 


(Pi  ~ A) 
(Pi  ~ A) 


( Pi  + A )( P/  - di)~2sr 


f (P,+  A )(P,~d,)  2 ^ 

5 3 5 


(P,~A)~  ^(P,-  A)  = P,  ~ A 


(-(P, -A)  + <0=  (P,  "A) 


(P,-d,)^ 
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Thus 


(A.3.8) 


-rt  (f,2)  = 


+ 


1 P^P-d.) 


vT 


UU1 


>i  y 


(Pi  ~st) 


( P,  + T X P,  -dj  2 
s 


\\ 


5; 


All1 


• yy 


+ 


r 


(Pi -Si) 


1 


(Pi-d,-) 


/, 


yy 


which  reproduces  the  result  (3.2.10)  for  the  interior  portion  of  the  directional  error  function.  For 
the  exterior  portion  g , = ( p l — d t)~  + {q  t — R)~  , we  express 

(A.3.9)  |r,(g,2)  = X-r,{p-d,f  +±r,(q,-R)2. 

We  note 

%niS\p-d,)  = uT, 

v (p,-d,)2  =2(p,~  d, ) V ',,;i  (i>l-d,)  = Hpl-dl)UT . 

Thus  by  the  Product  Rule, 

^ rt  (P,-d,  )2  = V ,v>7  | v ^ ( Pi  - </,  )2  = V XYZ  (( Pi  - di  )U 7 ) 

= (VxYz(Pi~di))UT  + ( P,  ~ di )Vxrzt/7 

so  that  by  (3.2. 1 ) and  (A.  1 .3), 


(A.3.10) 


~r,(Pl 


df  =A,Ut  +(Pi-d,)I 


Similarly,  we  derive  using  (A. 2. 1 ) 


Vl„Jqi-R)  = -^-Ur 

(9,  - R)2  = 2(9,.  - . (q,  -R)  = 2(9,  - /?) 


y x 

= 2 

r*-.l 

</.  2 

U J 

PV 
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We  proceed,  using  the  Product  Rule  as  well  as  (3.2. 1 ) and  (A.  1 .3): 


-r,(q,-R)2  = vm-v^(9l-R)2  = v 


XYZ 


( R 

\ 

— -1 

P,UT 

A Pi  y 

) 

f 

v 

v XYZ 

\ 

p,uT  + 

(y  xyz  Pi)uT  + 

V 

u J 

/ 

U-  J 

U J 

= R 


= R 


XYZ 


\ 

+ 

y xyz  p,)u'  + 

f*-ll 

U, 

) 

y 

u J 

— 0 V-PA ) 

v q\ 


\ 

/? 

> 

' /?  ) 

p,t/r  + 

1 

4i/r  + 

— 1 

/ 

J 

U J 

py  xyzu' 
py  xyzu' 

Pi*- 


Collecting  by  matrices, 


±r,(qi-R)2  = -?Z-uuT  +^4-a,u' 

2 4,  <h 


4.1/' r + 

--'1 

U J 

‘1,  J 

P,I 


RPi 

y 


f 2 


UUT  + 


/? 

\ 

( R ) 

H 

1 

a,ut  + 

q, 

) 

U J 

Pj ' 


and  substituting  (A. 3. 10)  and  the  above  into  (A. 3. 9)  gives 


\r,(.gf)=  -^lvut  + 

2 P, 


f 


R 

\ 

( 

R 

\ 

+ — 

-l 

AiUt  + 

P,  ~d,+ 

1 

Pi 

q, 

) 

\ 

U ) 

) 

Since 


, Rp 2 R 

1 + — C^  + — 


Rpl+R_  _ R(p;+q?)  _ R(X2  + Y2+Z 2) 


q,  q,  q,  q, 

we  have  finally 


q, 


q, 


(A.3.11)  X~r,{g;)=  -^UU‘  + 

2 q; 


T . RiX_l+y__+Zl)Aur  + 


q, 


Rp, 


\ Hi 


-d 


This  reproduces  the  result  (3.2.1  1)  for  the  exterior  error  function. 
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A.4  Gradients  and  Hessians  of  the  Orthogonal  Error  Function 

We  repeat  the  definitions  of  Chapter  4.  The  error  function  for  orthogonal  fitting  of  a sphere  is 
given  as 

Earth  = Yjhi  =2>,-*r  =wf  - 2Rwi  + R 

1=1  i=i 

where 


hi  = wi  — R 


wf  = (X-xi)2+(Y-yl)2+(Z-zi)2=(X2+Y2+Z2)-  2 d,p,  + d2 

with  the  notations  (3.1.1) 

Xi  = d,  cos  (pt  cos  9,  = d& , V,  = dt  sin  (p  t cos  Gi  = d Tji  , z.  = di  sin  <9,  = </,.£ 
and  the  definition  (3.1.2)  of  the  auxiliary  quantity  pt . The  vector 


~X-x~ 

~x~ 

p 

w = 

l 

y-y. 

- 

Y 

~d, 

n, 

Z-z,_ 

Z 

Jdi  _ 

= U - d.A, 


will  play  the  key  role.  Indeed,  all  of  the  following  gradients  are  multiples  of  W 


(A.4. 1)  Vx^w,2=2W, 


Vyizh;,  = VxizV^  = \ — ' W,  ’ 

2 VV  VV 


V 


xrz 


VC. 


— - — V w = - — W 

2 V XKZ  VV;  V rr  / 

w vv; 


Thus 


(A.4. 2) 


V h = V vv  = — W 

V XYZW  V XKZ  / I 

VV 


and  in  view  of 
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it  follows  that 


VxvA  = Vxvz(w, 


R)2  = 2 (w,  - fl)V 


XVZ 


w 


= 2W‘  — R\V 

I 

w 


(A.4.3) 


-V  h = 

~ v XKZ'S 


/? 


W 


v w/y 


Te  Hessians  of  the  quantities  considered  below  are  linear  combinations  of  the  two  matrices 


/ = identity,  and 


X-x, 

y - y, 
z-z, 


[X-x,  Y-y,  Z-z,]  = W,W,r 


In  particular, 

(A. 4.4)  HXYZw2  = Hxyz(X2  + Y2  +Z2)-2dHXYzPi  + H XYZd2  = H XYZ(X2 + Y2  + Z2)  = 21 
and  by  (A.  1 .2), 

--*-2WWr2  + - — 21 
4 w3  ' ' 2 vv. 

i 1 

As  ht  - wv  — R and  h 2 = w~  - 2 Rwt  + R2 , we  have 


T T 1 1 \ — 7 2 T” 7 T 2 1 1 y T 2 

n vv^W.  = TV  v.^w.-V  VV,W.  H ti  = — - 


Xre  rK/ 


4 w 


3 XYZ  rvi  XYZ 


Z W 


XKZ  ' i 


= — -wvv  + — / 

3 ' ' 

VV’  VC 


H XYZ^1,  H XYZWi 


— l-jwyv?  +—i. 

vv  vv. 

I I 


and 


'ID  'i  n 

HyY7h 2 = H yy7w2  -2RH  yy7w  = 21  + —WWr — -I 

XYZ.  i XYZ  i XYZ  i 3 1 1 


vc 


VV- 


so  that 


(A.4.5) 


" X>A  = 


f _R_ 

V w.j 


I+Xww’ 

vv: 
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A.5  Mixed  Derivatives  of  the  Orthogonal  Error  Function 

First,  we  differentiate  with  respect  to  the  range  variable  d,  , starting  with  the  key  quantity  w, 


(A.5. 1) 


Thus 


and 


3w2 

a 

= (X-  - 

dd, 

dd, 

dw, 

1 1 aw; 

ddt 

2 w,  aa, 

a i 

i aw 

aa  w, 

w;  aa 

(A  + Y “ + Z " — 2 dj  p,  + d ~ ) — — 2(  pi  — d, ) 


(~2)(Pi  -di)  = -—(pi  -d^ 

l W;  W 


j(-l )(Pi  -di)  = —{Pi  -d^ 

w w 


l l 


K = 3(w  _ Rf  = 2(W(  _ = _2^(pi  _dj)  = _2 

ad,  ad,  ad  w. 


V H';  J 


( Pi  ~ dj ) 


ly 

2 ^ aa 


2 / 


f _ P 

f /?  1 

-Mxyz- 

I 

"<3~ 

1 

rvT 

1-  — 

V wiJ 

^ XYZ  Pi 


so  that,  finally,  by  ( A.4. 1 ) and  (3.2. 1 ), 


(A.5.2) 


ly  })h: 


XYZ 


dd 


A lP,-d,W, 

W- 


r R A 


V WU 


A. 


We  move  to  differentiation  with  respect  to  the  bearing  variables  (p,.0r  We  apply  the  differential 
matrix  operator  (compare  Section  (A.3)) 

r,  =vm(v^J, 

to  the  individual  orthogonal  error  square  h~  = w2  - 2Rw,  + R~ , so  that 
(A.5. 3)  r^h;  ) = r,(w2)  -2fl/2(w,.). 

As  vi’2  = X2  + T2  + Z2  - 2a,./?,  + a2 , we  find  in  view  of  p,  = X£,  + Ft/,  + Z<^,  (see  (A. 3. 4)), 
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VL,?,  w.2  =-2d^ln,i,Pi  =~2d,U‘ 


w,  = = TT=WV5.,.£,  w?  = ~(r^.VT)  = -^W 


2/w f ""'s'  ' 2w, 


w 


Then  by  (A.4.1)  and  the  Product  Rule, 


0(«-r)  = Vxrz(vr(iniiwf)=Vxrz{-2diUT)=-2dVxrzUT  = -Id, l 


rM)  = VXK(V{,  ,,£,«'/)=  V 


XKZ 


y , x 
.ii,' 
v WI  y 


-^■VXKZ 


— U‘ 

v vv<  y 


= -</. 


\ 


V — 

v XYZ 


W. 


«y 


t/r  - —V  t/7'  = -<y 

17  v XYZ Ui 

W 


v w,  y 


£/7  -d—I 

w 


= d—WUT  -d—I , 
w w 


so  that  by  (A. 5. 3) 


r,(/i,!)  = -24/  -4  + 4— i 

w w 


-2d. 


/ /?  ' 


v vv<  y 


I -2d, -^WUT 
w 


And  finally 


(A. 5.4)  - r,(h~)  = -d , 


( R 


V W.J 


I -d,-^r\VUT  = -r/. 


/ 3 
W 


( ‘ R\ 


\\  wt  y 


R 

I +^WU‘ 


w 


' y 


Taking  into  account  again  the  orthogonality  relations  (3.2.10),  this  matrix  can  be  symmetrized  by 
substituting  W7=UT+Aj  for  U 1 in  the  above  expression,  yielding  the  matrix 


(A.5.5)  '-r,(h;)  = -d, 


R 


V WiJ 


I - d,—jWWT  = -d, 
w 


n R \ 


w 


w 


R 


I + —WW 


J T 


W 


33 


Acknowledgement 

We  are  grateful  to  our  NIST  colleagues  James  Filliben,  Alan  Heckert,  Joan  Rosenblatt,  and 
Nien-Fan  Zhang  for  information,  discussions  and  material  about  the  subject  of  variance 
propagation. 


REFERENCES 

[1]  W.  C.  Stone,  M.  Juberts,  N.  Dagalakis,  J.  Stone,  and  J.  Gorman,  "Performance  Analysis 
of  Next-Generation  LADAR  for  Manufacturing,  Construction,  and  Mobility,"  NISTIR 
vol.  7117,  2004. 

[2]  D.  E.  Gilsinn,  C.  Witzgall,  G.  S.  Cheok,  and  A.  Lytle,  "Construction  Object  Identification 
from  LADAR  Scans:  An  Experimental  Study  Using  I-Beams,"  NISTIR,  vol.  7286,  2005. 

[3J  P.  J.  Besl  and  N.  D.  McKay,  "A  method  for  registration  of  3-D  shapes,"  IEEE  Trans. 
Pattern  Analysis  and  Machine  Intelligence,  vol.  14,  pp.  239-256,  1992. 

[4]  R.  O.  Duda  and  P.  E.  Hart,  "Use  of  Hough  transformation  to  detect  lines  and  curve  in 
pictures,"  Comm.  ACM,  vol.  15,  pp.  11  - 15,  1972. 

[5]  K.  Kanatani,  Statistical  Optimization  for  Geometric  Computation.  New  York:  Elsevier, 
1996. 

[6]  A.  C.  Davidson,  Statistical  Models.  Cambridge:  Cambridge  University  Press,  2003. 

[7]  D.  A.  Freedman,  Statistical  Models,  Theory  and  Practice , revised  ed.  Cambridge: 
Cambridge  University  Press,  2009. 

[8]  P.  T.  Boggs,  J.  R.  Donaldson,  R.  B.  Schnabel,  and  C.  H.  Spiegelman,  "A  computational 
examination  of  orthogonal  distance  regression,"  Journal  of  Econometrics,  vol.  38,  pp. 
169  - 201,1988. 

[9]  P.  T.  Boggs,  J.  R.  Donaldson,  R.  H.  Byrd,  and  R.  B.  Schnabel,  "ALGORITHM  676 
ODRPACK:  Software  for  Weighted  Orthogonal  Distance  Regression,"  ACM 
Transactions  on  Mathematical  Software,  vol.  15,  pp.  348  - 364,  1989. 

[10]  W.  Gander,  G.  H.  Golub,  and  R.  Strebel,  "Least-squares  fitting  of  circles  and  ellipses," 
BIT,  vol.  34,  pp.  558-578,  1994. 

[11]  C.  M.  Shakarji,  "Least-Squares  Fitting  Algorithms  of  the  NIST  Algorithm  Testing 
System,"  Journal  of  Research  of  NIST,  vol.  103,  pp.  633-640,  1998. 

[12]  S.  J.  Ahn,  W.  Rauh,  and  H.  J.  Wamecke,  "Least-squares  orthogonal  distances  fitting  of 
circle,  sphere,  ellipse,  hyperbola,  and  parabola,"  Pattern  Recognition,  vol.  34,  pp.  2283- 
2303,2001. 

[13]  I.  D.  Coope,  "Circle  Fitting  by  Linear  and  Nonlinear  Least  Squares,"  J.  Optimization 
Theory  Applicat.,  vol.  76,  pp.  381-388,  1993. 

[14]  B.  A.  Jones  and  R.  B.  Schnabel,  "A  comparison  of  two  sphere  fitting  methods.,"  in  IEEE 
Proceedings  of  the  Instrumentation  and  Measurement  Technology , 1986,  pp.  104  -109. 

[15]  I.  Kasa,  "A  Circle  Fitting  Procedure  and  Its  Error  Analysis,"  IEEE  Trans.  Instrum.  Meas., 
vol.  25,  pp.  8-14,  1976. 


34 


[16]  Y.  Nievergelt,  "Computing  circles  and  spheres  of  arithmetic  least  squares,"  Comput. 
Phys.  Commun.,  vol.  81,  pp.  343-350,  1994. 

[17]  M.  Renault,  "Fitting  Circles  and  Ellipses  to  Data  Using  the  Least-Squares  Method,"  in 
http.V/www.  math. temple,  edit/- renault/ellipses,  html. 

[18]  H.  Spath,  "Least-Square  Litting  with  Spheres,"  J.  Optimization  Theory  Applicat.,  vol.  96, 
pp.  191-199,  1998. 

[19]  D.  Umbach  and  K.  N.  Jones,  "A  Few  Methods  for  Fitting  Circles  to  Data,"  IEEE  Trans. 
Instrum.  Meas.  , vol.  52,  pp.  1881-1885,  2003. 

[20]  C.  Witzgall,  G.  S.  Cheok,  and  A.  J.  Kearsley,  "Recovering  Circles  and  Spheres  from 

Point  Data,"  in  Perspectives  In  Operations  Research , F.  B.  Alt,  M.  C.  Fu,  and  B.  L. 

Golden,  Eds.  New  York:  Springer,  2006,  pp.  393-413. 

[21]  M.  Lranaszek,  G.  S.  Cheok,  K.  S.  Saidi,  and  C.  Witzgall,  "Litting  Spheres  to  Range  Data 
Lrom  3-D  Imaging  Systems,"  IEEE  Trans.  Instrum.  Meas.,  vol.  58,  pp.  3544-3553,  2009. 

[22]  M.  Franaszek,  G.  S.  Cheok,  K.  S.  Saidi,  and  C.  Witzgall,  "Sphere  Fitting  as  a Gauge  of 
the  Repeatability  of  3D  Imaging  Systems,"  in  preparation,  2009. 

[23]  GUM,  "Evaluation  of  measurement  data  - Guide  to  the  expression  of  uncertainty  in 
measurement  (http://www.bipm.org/en/publications/guides/gum.html ),"  JCGM  2008. 

[24]  H.  Ku,  "Notes  on  the  Use  of  Propagation  of  Error  Formulas,"  Journal  of  Research 
National  Bureau  of  Standards,  vol.  70C,  pp.  331  - 341,  1966. 

[25]  D.  M.  Himmelblau,  Process  Analysis  by  Statistical  Methods.  New  York:  Wiley,  1970. 

[26]  N.  R.  Draper  and  H.  Smith,  Applied  Regression  Analysis , 2nd  ed.  New  York:  Wiley, 

1981. 

[27]  C.  G.  Broyden,  "The  Convergence  of  a Class  of  Double-rank  Minimization  Algorithms: 
2.  The  New  Algorithm  " IMA  J.  Appl.  Math.,  vol.  6,  pp.  222  - 231,  1970. 


35 


